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Abstract 

The near-field screech-tone noise of a typical underexpanded circular jet issuing from a sonic nozzle 
is simulated numerically. The self- sustained feedback loop is automatically established in the simulation. 

The computed shock-cell structure, acoustic wave length, screech tone frequencies, and sound pressure 
levels in the near field are in good agreement with existing experimental results. 

1 Introduction 

Jet noise is an important and still challenging topic in aeroacoustics. Many of its aspects are of primary practical 
importance and the associated complicated physical phenomena are the topic of many experimental and theoretical 
investigations. References [1-4] provide a comprehensive discussion and further references. Under/over-expanded 
supersonic jets emit mixing noise, broadband shock-associated noise, as well as screech tones under certain conditions. 
The mixing noise is directly associated with large-scale structures, or instability waves, in the jet shear layer and the 
broadband shock-associated noise is caused by the interaction of these waves with the shock-cell structure in the 
jet core. The screech tones arise due to a feedback loop involving the large-scale structures developing in the jet 
shear layer, their interaction with the jet-core shock-cell structure producing upstream propagating acoustic waves, 
and re-generation of the large-scale structures at, or in the vicinity of, the nozzle lip. This feedback loop leading 
to screech tones is sensitive to small changes in the system conditions, and the understanding of the phenomena to 
date is based mostly on experimental observations [5-9] . Screech is of particular interest not only because of general 
noise-reduction concerns, but also because of potentially destructive structural interaction leading to sonic fatigue. 

Reliable numerical simulation of jet screech noise, i.e., near-field noise computation in the presence of shock cells 
in the jet core, has only quite recently become feasible. It is a challenge in that the numerical scheme must: (i) handle 
shock waves, (ii) resolve acoustic waves with low dispersion and dissipation errors, (iii) resolve the instability waves 
in the jet shear layer and their interaction with the core shock-cell structure, and (iv) have an effective non-reflecting 
boundary condition. Successful direct computation of screech for circular jets have been carried out by Shen and Tam 
[10, 12] using the Dispersion Relation Preserving, or DRP, scheme and by the present authors [13-16] using a Navier- 
Stokes solver based on the space-time conservation element and solution element, or CE/SE, method [17, 18]. In their 
3-D computation, Shen and Tam [12] adopted a spectral method in the azimuthal direction and by using only a limited 
number of spectral functions achieved substantial savings of computer memory and CPU time, without deterioration 
of accuracy. Other recent computational work includes [19, 20]. 

In general, the CE/SE method systematically solves a set of discretized space-time integral equations derived 
directly from the physical conservation laws and naturally captures shocks and other discontinuities in the flow. The 
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method may be categorized as a finite- volume method, where the conservation element ( CE ) is equivalent to a finite 
control volume (cell) and the solution element ( SE ) can be understood as the space-time cell interface. The scheme 
is second-order accurate in both space and time and the space-time surface fluxes are carefully calculated. Despite its 
formal second order accuracy, it possesses low dispersion error and low dissipation. Multi-dimensional non-reflecting 
boundary conditions (NRBCs), based on plane- wave propagation principles [21], are easily implemented. In principle, 
these NRBCs are equivalent to using characteristics but are simpler to apply. With these NRBCs, a small near field 
computational domain can be used. The CE/SE procedure is also easily adapted to complicated geometries using 
unstructured grids. The method is robust enough to cover a wide spectrum of compressible flow — from weak linear 
acoustic waves to shocks — hence, it is an appropriate tool for the present jet screech computation. Details of the 
method are described in [13, 17, 22]. 

2 The Axisymmetric Jet Screech Problem 

Figure 1(a) shows the geometry of the convergent nozzle in Panda’s experiment [9]. The flow is choked (i.e., Mach 
number M e = 1) at the nozzle exit. The operating condition of the jet is fully described by additionally specifying the 
plenum (reservoir) to ambient pressure and temperature ratios, respectively. The pressure ratio is commonly expressed 
in terms of the jet Mach number, Mj , which would be the exit Mach number of the corresponding perfectly expanded 
jet. For the cold-flow situation (unity temperature ratio) considered in Ref. [9], the dominating screech tone mode is 
axisymmetric for Mj ^ 1.19. Hence, the numerical simulation can employ an axisymmetric CE/SE Navier-Stokes 
solver for these operational conditions. Furthermore, the attention is focussed here on the near field of the nozzle since 
this is the noise source region. The diameter, D , of the jet nozzle is chosen as the length scale; the density, p^, 


14.3D 



Figure 1 . (a) Geometry of the convergent nozzle and flange in Panda’s experiment [9] (dimensions in mm), 
(b) Computational domain 


speed of sound, a oo, and temperature in the ambient flow are taken as scales for the dependent variables; and all 
other scales are constructed from these. The computational domain, see Fig. 1(b), spans between —8.3 < x ^ 6 and 
0 ^ r < 11.7, with x and r being the nondimensional stream wise and radial coordinates. The flow inside the nozzle 
is not computed, rather the steady flow conditions are prescribed at the nozzle exit which is located at x = 0. This 
inflow plane is recessed by two cells so as not to numerically restrict or influence the feed-back loop. The numerical 
scheme [13, 14] utilizes a triangulated grid of about 230,000 cells. The triangular cells are generated by dividing 
rectangle-like quadrilaterals. Non-reflecting boundary conditions are applied to the upper and outflow boundaries and 
a symmetry condition is applied at the center axis. The last 10 stream wise cells are exponentially stretched and serve 
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Figure 2. Numerical schlieren results (density gradient modulus) at t = 410,000 steps; Mj = 1.11, 1.15, and 1.19 in top, middle, 
and bottom rows; time-averaged and instantaneous results in left and right columns 


as a buffer, or sponge, zone to essentially eliminate any small remaining numerical reflection from the downstream 
outflow boundary. 

Initially, the entire flow field is at rest and at ambient conditions, i.e., (using nondimensional variables) 

Po = 1, Po = ~, To = 1, uo = 0, vo = 0, 

7 

where p, p , T, u, v, and 7 denotes the density, static pressure, temperature, streamwise velocity component, radial 
velocity component, and constant specific heat ratio, respectively. The jet flow is then impulsively started. At the inlet 
boundary, the conservative flow variables (see Appendix Eq. 2) and their spatial derivatives are specified to be those of 
the ambient flow, except at the nozzle exit, where an elevated pressure is imposed, i.e., the jet is under-expanded, as in 
the physical experiments. By using the ideal gas isentropic relations, it follows that the nondimensional flow variables 


NAS A/TM— 2003-2 1 2626 


3 




Figure 3. Experimental schlieren picture showing the shock-cell structure for Mj = 1.19 [9, Fig. 5(b)] 



Figure 4. Instantaneous differential isobars, at t=41 0,000 steps, showing screech waves; Mj = 1.11, 1, 15, and 1.19 in left, middle, 
and right panels 


at the nozzle exit, with M e = 1, are given by 


_ JPe 
Pe , 



2 + (7-l)M? 
7 + 1 


7 

7 — 1 


T, = 


2T r 
7 + 1’ 


= T}/\ 


V e = 0, 


where T r is the reservoir (plenum) temperature. We will also follow the experimental cold-flow condition where the 
reservoir temperature equals the ambient one, i.e., T r = 1. The Reynolds number Re = Da^/v = 570, 000, where 
v is the kinematic viscosity at the ambient conditions. The weighted average a-e CE/SE scheme [22] is used with the 
weighting parameter a = 1 and e = 0.5 in the present computation. 


3 Numerical Results and Comparison to Experiments 

In this section, various computed results are compared to the experimental findings of Panda [9]. Figure 2 displays 
numerical schlieren (density gradient modulus) contours, well after the start-up transients have passed out of the 
computational domain for the cases of Mj = 1.11, 1.15, and 1.19 — top, middle, and bottom rows, respectively. The 
left and right columns show time-averaged and instantaneous results, respectively. The figure clearly displays the 
shock-cell structure as well as its deformation and destruction by the spatially growing large-scale structures in the 
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jet shear layer. The shock-cell spacing of 0 .8D at Mj = 1.19 agrees well with the experimental result of Panda [9, 
Fig. 5(b)], reproduced here as Fig. 3. 

Figure 4 displays instantaneous numerical isobars well after the start-up transients have passed out of the com- 
putational domain for the cases of Mj = 1.11, 1.15, and 1.19. The contours shown are for the fluctuating part of 
the pressure, i.e., the time-averaged pressure at each spatial location has been subtracted. Distinct screech waves are 
observed to emit from the 3rd to the 5th shock-cell and are reflected at the flange/nozzle body. The screech wavelength 
(1 .6D) at Mj = 1.19 agrees well with the experiment [9]. Spectral analysis yields a computed screech frequency of 
8513 Hz, which agrees well with the experimental value of 8525 Hz. Figure 5 displays computed sound pressure 




(c) frequency (Hz) 


Figure 5. Narrow-band SPL distribution at nozzle exit (x = 0, r = 0.6), 75 Hz digital binwidth; (a) Mj = 1.11, Ai\ 8289 Hz; (b) 
Mj = 1 . 15 , Ai. 7542 Hz, A 2 : 9110 Hz; (c) Mj = 1 . 19 , A 2 : 8513 Hz 


level, SPL, distributions (narrow-band) at the nozzle exit lip wall (x = 0, r = 0.6) for Mj = 1.11, 1.15, and 1.19. The 
SPL distributions show distinct spikes corresponding to the A\ and A 2 axisymmetric screech-tone modes observed 
in Panda’s [9] experiment. For Mj = 1.15, both the A\ and A 2 modes are dominating features, whereas only the 
Ai mode is dominating at Mj = 1.11 and only the A 2 mode is dominating at Mj = 1.19. Note that in our earlier 
work [13], where a much simplified description of the nozzle external geometry was used, the A 1 -mode screech tone 
(not observed in [9]) was also prominent at the latter jet Mach number. This illustrates the sensitivity of the screech 
phenomenon to geometry changes. Figure 6 shows a comparison of the current numerical with Panda’s experimental 
results for the Ai and A 2 screech-tone frequencies and SPL at the nozzle lip, both as functions of the jet Mach number. 
The screech-tone frequencies are well predicted, whereas the numerical SPL values are somewhat low. 
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(a) M j (b) 


Figure 6. (a) Screech-tone frequency vs. Mj. (b) Screech-tone SPL at nozzle exit (x = 0 , r = 0 . 6 ) vs. Mj. Square: A\ mode, 
circle: A 2 mode; filled symbol: numerical simulation, open symbol: Panda’s experiment 


The SPL along an inclined line at the outer edge of the shear layer is shown in Fig. 7 for Mj = 1.19. This 
figure shows Panda’s [9] data for the A 2 axisymmetric screech mode, the corresponding result from the simulation, 
as well as the computed total SPL and subharmonic of the A 2 mode. Even though the A 2 SPL level in the vicinity 



Figure 7. Comparison of computed and experimental [9] SPL along the shear layer edge, Mj = 1.19 


of the nozzle lip is too low, its early stream wise growth rate is, however, well predicted indicating that the jet shear 
layer is in general well resolved. Nonlinearity in the shear layer limits the amplitude that can be obtained and once 
the SPL level ‘catches up’ the agreement is very good, in particular the SPL level in the stream wise region where 
the backward radiating acoustic waves are generated is well predicted. The strong subharmonic that appears further 
downstream is due to (axisymmetric) vortex pairing in the shear layer and as a consequence a second streamwise peak 
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at the A 2 frequency occurs due to nonlinear effects. However, 3-D effects are most likely to have come into play in the 
experiment at these streamwise locations leading to suppression of the subharmonic. The subgrid-scale eddy-viscosity 
model used in the present computation (see Appendix Eq. 3) does reduce the subharmonic somewhat, but cannot 
supress it since the model does not describe the effects of large-scale 3-D motions that ultimately become important 
at these streamwise locations. The reason for the low SPL values at the nozzle lip in the simulation is currently being 
pursued. 

Figure 8 shows computed SPL contours for the A 2 mode. As in the experiment [9], a standing wave structure 
can be observed along the edge of the jet shear layer. As explained by Panda [9], the standing wave formation is due 
to combination of the downstream-propagating hydrodynamic and upstream-propagating acoustic fluctuations, both 
integral parts of the screech feedback loop. The results are in general agreement with the experimental ones except, 
as pointed out earlier, that the computed SPL levels are too low in the vicinity of the nozzle lip and the existence of a 
second, or extended, region of elevated values further downstream ( x/D > 4). 
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Figure 8. Computed A 2 -mode SPL levels for Panda’s experiment [9], Mj = 1.19 


4 Concluding Remarks 

In this paper, screech-tone noise analysis of an underexpanded supersonic axisymmetric jet was carried out using 
a finite- volume numerical scheme. Many hydrodynamic and acoustic aspects of the computed results are in good 
agreement with experimental findings [8, 9], including the variation of screech frequency with Mach number for the 
axisymmetric Ai and A 2 modes. Since no forcing was applied in these simulations, the numerical results clearly 
indicate the presence of a self- sustained oscillation. LES modeling was used to better capture the jet spreading further 
downstream, i.e., after the initial shear-layer roll up, but is not believed to be a critical element for the screech feedback 
loop. Note that screech is essentially an inviscid phenomenon. It is concluded that the simulation shows a reasonably 
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good agreement with experimental data in the streamwise region where the flow is expected to be predominantly 
axisymmetric and, hence, that jet screech is successfully simulated using the present scheme. 

The advantages of the present scheme were confirmed in the simulation. The implementation did not require any 
special treatment or parameter selections to capture both the shock-cell structure nor the acoustics. The robust finite- 
volume NRBC in combination with buffer zones virtually eliminates numerical reflections and allows for an effective 
near-field computation. The use of an unstructured grid also offers flexibility for dealing with more complicated 
geometries. 
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Appendix: Governing Equations 

The axisymmetric Navier- Stokes equations can be written in the following vector conservation form: 

Ut + F x + G r = Q, (1) 

where x, r > 0, and t are the streamwise and radial coordinates and time, respectively. The four components of the 
conservative flow variable vector U are given by 

Ui=p, U 2 = pu, U 3 = pv, U 4 =p/('f- 1) + p(u 2 + v 2 )/2. (2) 

The components of the flux vectors F and G are given by 

F\ = U 2 , F 2 = ( 7 - 1 )U 4 + [(3 - 7 )Ul - (7 - l)Uf\ /2U 4 - p(2u x - §A), F 3 = U 2 U 3 /U 4 - p(u r + v x ), 

'"y c) Ua r ii^‘ 

F 4 = - (7 - 1 )U 2 [Ui + Ui] /2 Ui - p 2 uu x + v(u r + v x ) - fuA + |— ) , 

G ! = U 3 , G 2 = U 2 U 3 /U\ - p(u r + v x ), G 3 = ( 7 - 1 )U 4 + [(3 - 7 )C/ 3 2 - (7 - 1 )U%] /2 U x - p(2v r - §A), 

ry T J A 7/^ —I— 7 ;^ 

G 4 = 7 U 3 u 4 /u 4 - (7 - 1 )U 3 [Ui + Ui] /2 Ui - p 2vv r + u(u r + v x ) - frA + |— ) , 

where u x ,u r ,v x , v r are the spatial derivatives of the streamwise and radial velocity components, which can be written 
in terms of the conservative variables C/2, Us and C/4, with Pr (= 0.72) being the Prandtl number, fi = 1/ Re the 
nondimensional viscosity, and the velocity divergence 

A = u x + v r + v/r. 

The four components of the ‘source’ term Q are given by 

Qi = ~U 3 /r, Q 2 = -U 2 U 3 /U ir , Q 3 = -Ui/U 4 r, Q 4 = -G 4 (p = 0)/r. 

The above nondimensional equations would form the basis for a direct numerical simulation (DNS). Currently, 
this is not fully feasible for the high Reynolds number flows of interest here. The present computations are of the 
large-eddy- simulation (LES) type [23], where scales smaller than the grid resolution are modeled. In this case, the 
equations above are considered the filtered ones governing the resolved scales, where p and p are to be interpreted as 
simply averaged over the subgrid scales and the other dependent variables are to be interpreted as Favre, i.e. density- 
weighted, averages. The simplest model for the effects of the unresolved scales on the motion is obtained by using a 
Boussinesq eddy-diffusivity assumption for the subgrid shear stresses and heat flux coupled by Smargorinsky’s model 
for the eddy viscosity (eg. [24]), i.e., p above is replaced by 

p = l/Re+(C s S) 2 (S ij S ij ) 1/2 , (3) 
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where Sij = \ ^ is the rate-of-deformation tensor, 5 = (dxdr) 1 / 2 is a local measure of the grid size, and 

C s is a (nondimensional) constant. Note that the simplifying assumption that the subgrid-scale Prandtl number equals 
the laminar one is also made here — if it is not, then the G 4 component needs to be slightly modified. 

Panda et al. [25] investigated the initial status of the shear layer at the nozzle exit using a molecular Rayleigh 
scattering technique for a range of Mach numbers that spans the current one of interest and found it to be nominally 
laminar. This means that the small-scale turbulence that needs to be modeled is mainly produced as the jet shear layer 
rolls up, hence it is limited to a certain region of the spatial computational domain. Consequently, in the present work, 
C s is essentially set to zero for r > 1.5 or x < 3 and is linearly ramped up for 3 < x < 5 to the value C s = 0.1 
(see [23]), which is then used for r < 1.5 and x > 5. 

By considering (x, r, t) as coordinates of a three-dimensional Euclidean space, E%, and using Gauss’ divergence 
theorem, it follows that Eq. (1) is equivalent to the following integral conservation law: 

(f H m -dS = f Q m dV , m = l,2,3,4, (4) 

JS(V) Jv 

where S(V) denotes the surface around a volume V in E% and H m = (F m , G m , U m ). These equations are discretized 
and solved using the a — e weighted-average CE/SE scheme with a = 1 and e = 0.5, see [13, 14, 22] for further details. 
Note that Fureby and Grinstein [26] have pointed out that applying flux limiters to finite- volume methods, even in the 
absence of any explicit LES assumptions, effectively leads to LES schemes with minimal implicit SGS models. They 
demonstrated through ‘error’ analysis of a particular scheme that the flux-limiters (essentially low-pass frequency 
filters) build into the algorithm produces additional terms in the equivalent differential forms of the momentum and 
energy equations that can be interpreted as the SGS stress tensor and flux, respectively. Hence, there is an additional 
implicit SGS model inherent in the particular CE/SE scheme used in addition to the explicit assumption in Eq. 3. 

The computational domain was decomposed into 10 subdomains using the METIS mesh partitioning code [27] 
freely available from the University of Minnesota. The CE/SE solver implemented MPI [28], or message passing 
interface, calls to exchange pertinent data between neighboring subdomains. The parallel computations were carried 
out on the CW-7 Linux PC (PHI) cluster at NASA Glenn Research Center. 
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